Non-equilibrium phase transition in a two-temperature lattice gas 



(N 



T3 
O 



> 



0\ 



- 1—1 

X 



Attila Szolnoki 

Research Institute for Materials Science, H-1525 Budapest, POB 49, Hungary 

t> ' 

' A two-temperature lattice gas model with repulsive nearest-neighbour interactions is studied 

i using Monte Carlo simulations and dynamical mean-field approximation. The evolution of the two- 

■ dimensional, half-filled system is described by an anisotropic Kawasaki dynamics assuming that 

— ' the hopping of particles along the principal directions is governed by two heat baths at different 

temperatures T x and T y . The system undergoes an order-disorder phase transition as T x (T y ) is 
' varied for sufficiently low fixed T y (T x ). The non-equilibrium phase transition remains continuous 

and the critical behaviour belongs to the Ising universality class. The measure of violation of the 
fluctuation-dissipation theorem can be controlled by the value of the fixed temperature. We have 
found an exponential decay of spatial correlations above the critical region in contrast to the two- 
l ' temperature model with attractive interactions. 
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I. INTRODUCTION 

> \ 

, t \ The study of driven diffusive lattice gas model, which served a convenient tool to produce stationary non-equilibrium 
state, revealed a lot of new phenomena compared to the corresponding equilibrium properties OJ . Such a phenomenon 
is the emergence of a novel non-equilibrium fixed point of critical behaviour in the case of attractive interactions or 
the surprising feature of the fluctuations of internal energy [^]||] . In this model an external uniform electric field is 
applied to induce a preferential hopping in one direction leading to a permanent particle transport in the stationary 
state [Q. Unfortunately, it is rather difficult to realise the model in laboratory because the driving field cannot be 
q ■ derived as the gradient of an electrostatic potential. 

This problem can be avoided if we use an alternating field instead of the uniform one, which also enlarges the 
probability of jumps along specified lattice directions without resulting in a non-zero average particle current. Similar 
effect can be obtained in a two-temperature model where the hopping of particles along the principal directions is 
. governed by two heat baths at different temperatures |||J . Now the flow of energy between the heat baths drives 
the system into a non-equilibrium stationary state. The non-equilibrium phase diagrams have already been studied 
£N1 ■ in two-dimensions in the case of attractive interactions Q| and for nearest-, and next-nearest neighbour repulsive 
ON ■ interactions Common feature of the two diagrams is the existence of an ordered phase even at T x = oo, if T y is 
low enough. Similar behaviour is observed when the uniform electric field drives the systems into non-equilibrium 
states. Here the ordered phases also exist at infinite E if the temperature is low enough Completely different 

behaviour is found for nearest-neighbour repulsive interactions where the field induced interfacial transport prevents 
the long-range ordering JlO| . Although similar effect cannot be expected for the corresponding two-temperature model 
' because of the zero macroscopic current, still the stability of the ordered phase is doubtful Q|. 

In this paper we determine the phase diagram of a two dimensional two-temperature model with nearest-neighbour 
repulsive interactions. For this purpose we perform dynamical mean- field analysis to calculate the phase transition 
points. These results are supported by Monte Carlo simulations. We also study how the non-equilibrium conditions 
influence the critical behaviour of the system. This problem is especially exciting because similar non-equilibrium 
perturbation results in a critical behaviour different from its equilibrium counterpart in some models where the 
global current vanishes PHI- The finite-size scaling of the MC data indicates that the present model falls into the 
universality class of the equilibrium Ising model in agreement with a theoretical prediction |]l2| . 

The next section is devoted to the violation of fluctuation-dissipation theorem (FDT) which is generally one of the 
benchmarks of non-equilibrium phenomena. The breakdown of this theorem results in a difference between the specific 
heats defined as the temperature derivative of the average energy on the one hand and derived from the fluctuations of 
energy on the other hand. For example, in the uniformly driven model with attractive interaction the field suppresses 
the fluctuations and the difference is relevant [p]J^Jl^|. Using both ways we calculate the specific heat for our model 
and show that the rate of the violation of FDT can be influenced by the variation of the fixed temperature. 

A wide variety of non-equilibrium systems exhibit " generic scale invariance" , namely, long-range spatial correlations 
for arbitrary parameter values. An example of such a system is the two-temperature model with attractive interactions 
where a power law decay of pair correlation function is observed for all temperatures above criticality |J . We also 
examine the two-point correlations of our model and find short-range correlations at high temperatures. 
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II. THE MODEL 



Consider a two-dimensional lattice gas on a square lattice with L x L = N sites under periodic boundary conditions. 
The occupation variables n, = or 1 (if the lattice site is empty or occupied) and we assume half-filled occupation 
(Si n * = The Hamiltonian is given by 

H = -jJ2n i n j (1) 

where the summation is over nearest-neighbour pairs and J = — 1. Kawasaki exchange dynamics is assumed [ jl4) , 
namely the particles can jump to one of the empty nearest-neighbour sites. Anisotropic hopping rates are used for 
directions a = x or y: 

gaiAH) = 1 + cJaH/T.) (2) 

where AH is the energy difference between the final and initial configurations. The lattice constant and the Boltzmann 
constant are chosen to be unity as usual. 

In the case of T x = T y , the above model is equivalent to the celebrated kinetic Ising model which undergoes 
a continuous phase transition at T/v = 0.567. Below Tjv the system orders into a checkerboard-like pattern. If 
the lattice is divided into two interpenetrating sublattices (A and B) then particles prefer staying in one of these 
sublattices. Hence the difference between the two averaged sublattice occupations is the global order parameter: 
m = \pa ~ Pb\i where pa denotes the averaged occupation of sublattice A. A suitable local order parameter is the 
staggered magnetization at site j: cf>j = (— iy (2rij — 1). It is important to note that <j) is not conserved resulting in 
an important effect on the spatial correlations at high temperatures. 

When T x and T y are different the energy flows permanently from the heat bath of higher temperature to the 
lower one through the system and maintains a stationary non-equilibrium state. In other two-temperature model the 
previous studies set one of the temperatures (e.g. T y ) to be infinite, which reduced the number of control parameter 
to one ||[l5| and the phase transition was studied as a function of T x . Now we cannot use this simplification because 
the long-range ordered phase is expected to be unstable when one of the temperature is too high (the explanation of 
this expectation will be given in the next section). Therefore we keep T y at a constant (finite) value when we vary T x 
to find the order-disorder transition point. Naturally similar behaviour would be obtained by exchanging the roles of 
the two temperatures. 



III. PHASE TRANSITION AND CRITICAL BEHAVIOUR 



The study of non-equilibrium phase diagram by dynamical mean- field approach proved to be useful in many instances 
H[l6|. In this approach we determine the stationary solution of the equations describing the time evolution of 
configuration probabilities. At the level of /c-point approximation, these probabilities are given by probabilities of the 
possible configurations on fc-point clusters. We refer to [|l6) for further details of this method. 

First we use the two-point level, consequently correlations larger than two lattice constants are neglected. In spite 
of this simplification, the simple approximation predicts the relevant features of the phase diagram well as shown in 
Fig. |l|. Namely, the ordered phase becomes unstable even at T y = if T x is high enough. This behaviour may be 
understood qualitatively because the nearest-neighbours of the occupied lattice sites are empty in the ordered phase. 
Therefore the high temperature can influence the jumps and particles may leave the preferred sublattice. The phase 
transition remains continuous at any value of T y . Evidently, the quantitative predictions of this approximation are 
far from the correct values (e.g. T? P (T X = T y ) = 0.72, T? p (T y = 0) = 1.18). 

To obtain more accurate data, we apply four-point approximation. Here we have to solve a system of eight coupled 
nonlinear equations numerically. The results (e.g. T^ P (T X = T y ) — 0.60, T^ p (T y = 0) = 0.76) are already comparable 
to the data of Monte Carlo simulation as shown in Fig. [l| (The simulation data T^' IC (T X = T y ) = 0.56 and 
(T y = 0) = 0.74)). The phase transition is continuous even at T y — 0, therefore the order parameter decreases 
gradually when T x tends to the the transition point as Fig. [2] shows. Naturally, we have the same function when 
T x — and T y is varied. As a consequence the system becomes not perfectly ordered when we cool it "horizontally" 
to T x = axis on the phase diagram at a fixed nonzero value of T y . Namely, the order parameter saturates to a value 
less than 1 when T x — > 0. The visualisation of particle distributions during the MC simulation shows that there are 
still some defects in the ordered homogenous phase even at T x — because some particles do not stay in the preferred 
sublattice in consequence of the nonzero value of T y . 
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The MC simulations were performed on a square lattice (L x = L y = L) where system sizes ranged from L = 20 to 
150. At given (T x , T y ) temperatures following a thermalization of 2, 000 — 40, 000 MC steps, the time interval between 
two independent measurements ranged from 5 MCS for L = 20 systems to 100 MCS for L — 100 systems. Increasing 
T x the whole procedure was repeated. To preclude the possibility of first-order transition, especially at low T y , we 
repeated the measurement starting from high T x and decreasing T x , but we have not detected hysteresis in the order 
parameter function. It is worth mentioning that the anisotropic dynamics does not result in anisotropy, namely the 
pair correlation functions to the x and y directions are identical within the limit of accuracy. 

We turn now to study the critical behaviour of the model. For the previously studied two-temperature models 
the non-equilibrium behaviours differ from the equilibrium counterparts. In case of attractive interaction instead of 
the equilibrium Ising universality class, the anisotropic finite-size scaling resulted in new exponents (3 = 0.33 and 
v = 0.60 jig ]. The other model with repulsive nearest- and next-nearest neighbour interactions leaves the equilibrium 
universality class of the x-y model with cubic anisotropy and is described by the Ising exponents ||] . Common features 
of these models are that the degenerate ground states violate the x-y symmetry of the systems. For example in case 
of the attractive interaction the particles condense into a strip oriented either horizontally or vertically in the ordered 
phase. In the other previously studied model with repulsive interactions the particles form alternately occupied and 
empty columns (or rows). The fact that the anisotropic dynamics breaks the x-y symmetry and reduces the number 
of the stable groundstates may result in a new non-equilibrium behaviour. In contrast to these models, in the present 
model the checkerboard-like groundstate does not break the x-y symmetry, therefore the non-equilibrium perturbation 
does not reduce the number of the possible groundstates. 

As a consequence of the previous argument our conjecture is that the critical behaviour of the non-equilibrium model 
may be described by Ising exponents. This expectation agrees with the prediction of the field theoretic renormalisation 
group analysis 

Although the dynamical mean-field approach provides a correct phase diagram, it is not applicable to describe 
critical behaviour. To study the critical behaviour we have performed a finite size scaling of MC data where the usual 
scaling form of the order parameter was used: 

m{L,T x ) = L-P/ v m{T x L l / v ) . (3) 

Here t x denotes the reduced horizontal temperature deviation from the critical temperature. Notice that T y is 
constant. Using (3 = 1/8 and v = 1 in the scaling form, we have got nice data collapse at T y = as Fig. || shows. 
Similar excellent data collapses were found when we repeated the finite-size scaling at T y = 0.33 and 0.64. The Ising 
universality exponents are believed to hold along the entire line of phase transitions. 



IV. NON-EQUILIBRIUM FLUCTUATIONS 



The character of the phase diagram makes possible to study the phase transition in a wide range of T y . The non- 
equilibrium phase transitions are generally accompanied by the violation of the FDT. It results in difference between 
the specific heats defined as the derivative of energy with respect to the temperature (Cd = dE/dT) and derived 
from the fluctuations in the energy (Cf = ({E 2 ) — (E) 2 )/kT 2 , where ( ) stands for averaging over time). To check 
this Cd and Cp are compared at different fixed values of T y when we cross the phase transition line by increasing T x . 
In this case the correct definition of the two quantities arc: 

° D -dT x 

At T y = the fluctuations are suppressed by the low temperature, therefore Cf is smaller than Cd even at the 
phase transition point. Increasing T y , larger fluctuations can appear, which decrease the difference between Cd and 
Cf- At a special value of T y = 0.33, the specific heats, which are calculated different ways, collapse. Further increasing 
T y , the fluctuations become larger and Cf exceeds Cd- Figure ^indicates, how strongly the difference depends on T y . 
Thus we can control the rate of violation of FDT by changing T y , moreover, FDT seems to be valid at a special value 
of T y . It is easy to detect that the Cf curve is not symmetrical around the transition point (T c ) when T y = 0.64. This 
feature may be understood qualitatively. When T x is increased at fixed value of T y near 0.64, the system is close to 
the phase transition line in the sub-critical region (T x < T c ). While in the super-critical region (T x < T c ) the system 
is further away from the critical line. Therefore the fluctuations are enhanced in the T x < T c region comparing to the 
T x > T c region. 

Though the data are obtained from small systems, similar behaviour can be observed for larger systems. This fact 
is illustrated in Fig. || where the specific heats are plotted for different system sizes at the above mentioned special 



and Cf = 



m (e)< 

kT 2 



(4) 
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value of T y . It is worth mentioning that in the corresponding uniformly driven diffusive model the sharpening of 
the specific heat is not a monotonic function of the system size fl7| . The reason of this unusual behaviour is that 
the self-organising domain structure, which characterises the thermodynamic limit, can only appear when the system 
sizes is large enough. We expect traditional size dependence of specific heat in the two-temperature model because 
the reflection symmetric dynamics does not cause enhanced interfacial transport here. Although the peaks of the 
specific heats are more moderate than in the equilibrium case, the intensity of the peaks increases with lattice size 
suggesting a divergence in the thermodynamic limit as Fig. [s| shows. 

Finally, we consider the spatial correlations above the phase transition point. Long-range correlations are so often 
observed in non-equilibrium systems that the exponential decay of spatial correlations seems to be exception. It is 
believed that in the non-equilibrium steady-state the spatially anisotropic conserved dynamics produces such long- 
range correlations jl^,^| ■ In agreement with this argument both simulations |jj and perturbative approaches |19| ] found 
power law decay of correlations at high temperatures in different two-temperature models. In our model the order 
parameter, which is the 'staggered' density, obeys no conservation law therefore the theory predicts exponentially 
decaying correlations. However, the question comes up whether the validity or the violation of FDT may influence 
the spatial correlations somehow at high temperatures in the present model. 

To measure the equal time pair-correlation function 

G(r) = <#O)0(r)) , (5) 

where <f> is the local order parameter, we chose the system size to be considerably larger than the correlation length. 
Thus we have chosen L ranged from 50 to 100 when we measured G(r) at T y = 0.33 for three different T x which are 
above the transition point (T x = 1.2T C , 1.3T C , and 1.4T C ). In every case G(r) converged to zero. We have plotted 
log G(r) vs. r in Fig. to demonstrate the exponential decay of correlations. Obviously the negative slopes become 
smaller when T x decreases resulting in larger correlation lengths. Similar behaviours were found for other values of 

Ty. 



V. CONCLUSION 



The phase transition and critical behaviour of a two-temperature model with repulsive nearest-neighbour inter- 
actions were studied by dynamical mean-field approximation and Monte Carlo simulations. According to both 
approaches, the system undergoes a phase transition if one of the temperatures is high enough. The phase tran- 
sition is continuous even at zero fixed temperature and the critical behaviour is not affected by the non-equilibrium 
perturbation. This result agrees with the earlier theoretical expectations. 

The present model gives an example that FDT can be valid or violate in case of the same non-equilibrium model at 
different values of the control parameter. The Monte Carlo data of the pair correlation function resulted exponential 
decay in the high temperature regime, which is in agreement with the theoretical prediction since the order parameter 
is not conserved. We found that the nature of spatial correlations are independent of the violation FDT. 
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FIG. 1. Phase diagram predicted by two-point (dashed line) and four-point (solid line) dynamical mean-field approximations. 
Filled boxes indicate the results of Monte Carlo simulation. 
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FIG. 2. Order parameter as a function of T x suggested by four-point approximation at T y — . 



1 1 


1 1 r— 




ordered \ 






•5 




FIG. 3. Finite size scaling plots of MC data for the order parameter at T y — 0. System sizes are L — 20 (x), 30 (□), 50 (A), 
and 100 (•). The slopes of the inserted lines are 1/8 and —7/8. 




FIG. 4. MC data for specific heats at different values of T y (0.64 (x), 0.46 (A), 0.33 (□), and (•). Solid lines represent the 
specific heat calculated from the energy fluctuation. The system sizes are 20 x 20. 
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FIG. 5. Specific heats calculated both ways at T y = 0.33 for different system sizes L = 20 (•), 30 (o), 50 (+), and 100 ( □). 
The error bars of Co are indicated. 
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FIG. 6. Pair correlation functions at different T x above the critical region when T y — 0.33 for all cases. T x = 1.2T C (A), 
1.3T C (•), and 1.4T C (□). 
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